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Summary. Granular materials fiuidized by a rapidly vibrating bottom plate often 
develop a fascinating density inversion: a heavier layer of granulate supported by a 
lower-density region. We employ the Navier- Stokes granular hydrodynamics to follow 
a density inversion as it develops in time. Assuming a dilute low-Mach-number flow, 
we derive a reduced time-dependent model of the late stage of the dynamics. The 
model looks especially simple in the Lagrangian coordinates. The time-dependent 
solution describes the growth of a density peak at an intermediate height. A transient 
temperature minimum is predicted to develop in the region of the density peak. 
The temperature minimum disappears at later times, as the system approaches the 
steady state. At late times, the predictions of the low-Mach-number model are in 
good agreement with a numerical solution of the full hydrodynamic equations. At 
an early stage of the dynamics, pressure oscillations are predicted. 



1 Introduction 

An ensemble of hard spheres in motion, that collide inelastically and are 
characterized by a constant coefficient of normal restitution, represents the 
simplest model of granular gas [1]. It also provides an excellent example of 
complexity. Indeed, the properties of the granular g whole, such as clus- 

tering [2, 3, 4, 5], come out as non-trivial consequences of the simple algebraic 
laws that govern the velocities of the binary collisions of individual particles. 
Furthermore, being intrinsically far from thermal equilibrium, granular gas 
defies statistical mechanics, the cornerstone of the equilibrium science. As the 
result, formulating a universally applicable continuum theory of granular gas, 
even in the dilute limit, is not a simple task [6, 7] . The use of the Boltzmann (or 
Enskog) equation, the starting point of a systematic derivation of the consti- 
tutive relations, is based on the Molecular Chaos hypothesis. This hypothesis 
is justified, in the dilute limit, for an ensemble of elastic hard spheres. Its use 
for inelastic hard spheres is, rigorously speaking, an uncontrolled assumption. 
Indeed, inelasticity of the particle collisions introduces inter-particle correla- 
tions which can invalidate the Molecular Chaos hypothesis. The correlations 
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become more pronounced as the inelasticity of the collision increases. On the 
contrary, for nearly elastic collisions, 1 - r 2 < 1 (where r is the coefficient of 
normal restitution) the correlations are small, and the Boltzmann equation 
can be used safely (again, in the dilute limit). 

An additional crucial assumption, made in the process of the derivation 
of the hydrodynamics from the Boltzmann equation, is scale separation. The 
hydrodynamics is valid if the mean free path of the particles is much less than 
any characteristic length scale, and the mean time between two consecutive 
collisions is much less than any characteristic time scale, described hydrody- 
namically. This assumption (as well as the dilute gas assumption, if used) 
should be verified, in every specific system, after the hydrodynamic problem 
is solved and the characteristic length and time scales determined. One more 
criterion for the validity of granular hydrodynamics is the smallness of fluctu- 
ations, driven by the noise caused by the discrete nature of the particles. The 
accuracy of hydrodynamics improves with an increase of the total number of 
particles in the system N [8] . 

Once hydrodynamics is valid, it is highly rewarding using it. Hydrody- 
namics has a great predictive power and helps to identify important collective 
phenomena (shear flows and vortices, shocks, different modes of clustering 
flows, phase separation etc.) that are impossible to perceive in the language 
of individual particles. Once identified, these collective phenomena can then 
be investigated in experiment and simulations in more general (not necessar- 
ily hydrodynamic) formulations. A recent example: employing a variant of 
granular hydrodynamics, Livne et al. [9] predicted a novel phase separation 
instability in a system of inelastic hard spheres, driven by a rapidly vibrating 
plate at zero gravity. This system had been investigated earlier by many work- 
ers. However, the phase separation instability, and a plethora of interesting 
effects accompanying it, were overlooked. Undoubtedly, the reason for that 
was that hydrodynamics had not been fully appreciated and exploited [10]. 

In this paper we employ the Navier-Stokes granular hydrodynamics for 
a theoretical analysis of the process of formation of a density inversion in a 
dilute granular gas driven by a rapidly vibrating bottom plate in a gravity 
field. The phenomenon of a density inversion in the steady state of this sys- 
tem is well known in the granular community. On the contrary, the dynamics 
of the formation of the density inversion has never been investigated previ- 
ously. As we will see, this dynamics is quite instructive. The hydrodynamic 
theory yields predictions of time-dependent quantities that can be checked 
in molecular dynamics simulations and experiment. Our additional motiva- 
tion here is pedagogical. While dealing with this problem, we will present two 
useful hydrodynamic techniques. The first is reduction of the order of the hy- 
drodynamic equations, based on a hierarchy of the length/time scales in the 
problem. The second technique is the use of the Lagrangian coordinates. 

Here is a layout of the remainder of this paper. In Sec. 2 we will write down 
a set of hydrodynamic equations and boundary conditions that describe the 
dynamics of formation of a density inversion. By introducing scaled variables, 
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we will delineate the two scaled governing parameters of the problem. In Sec. 
3 we briefly review the properties of the steady state and the condition for 
a density inversion. In Sec. 4 we assume a low-Mach-number flow, derive a 
reduced hydrodynamic model and use it for an investigation of the late-time 
dynamics of the formation of the density inversion. Section 5 reports a numer- 
ical solution of the full hydrodynamic equations and compares the numerical 
results with the predictions of the low-Mach-number model. A discussion of 
the results is presented in Sec. 6. 

2 The model problem and hydrodynamic equations 

We will adopt the simplest possible model, in two dimensions, that exhibits 
a density inversion. Let N 1 identical nearly elastic hard disks of diameter 
d and mass m = 1 move without friction in a two-dimensional box of width 
L x and infinite height. The driving base is at y = 0, the side walls are at 
x = and x = L x . The gravity acceleration g is in the negative y direction. 
The rapid vibrations of the base can be modeled, in a simplified way, by 
prescribing a constant granular temperature T = T at y = (we measure the 
granular temperature in the units of the velocity squared). This means that, 
upon collision with the base, the particle velocity is drawn from a Maxwell's 
distribution with temperature T . The kinetic energy of the particles is being 
lost by inelastic hard-core collisions parameterized by a constant coefficient of 
normal restitution r. Particle collisions with the side walls are assumed elastic. 
Alternatively, one can specify periodic boundary conditions in the x-direction. 

Introducing a hydrodynamic description, we deal with three coarse-grained 
fields: the particle number density n(r, t), mean flow velocity v(r, t) and gran- 
ular temperature T(r, t). In addition, we make two important assumptions. 
First, we assume that that the collisions are nearly elastic: 1 — r 2 <C 1. Second, 
we assume that the granular gas is dilute everywhere (including the density 
peak region, see below), that is n is small enough compared to the close- 
packing density n c . These assumptions enable us to use the classic version 
of the Navier-Stokes granular hydrodynamics [11], and to limit ourselves to 
the leading order terms in n/n c in the constitutive relations. The dilute gas 
condition will be checked a posteriori. 

It has been shown recently that, when this system is large enough in the 
lateral direction, thermal convection may occur [12, 13, 14]. In this work we 
will assume that L x is sufficiently small, and the convection is suppressed by 
the viscosity and lateral heat conduction. Therefore, we assume flow solely in 
the vertical direction. Therefore, the governing equations are: 




(1) 




4 Y. Bromberg, E. Livne, and B. Meerson 



fdT ^dT\ _ dq p dv 
\dt dy J dy dy 

Here P = —nT + ndv/dy is the vertical component of the stress tensor, 
q = —ndT/dy is the heat flux, and I is the heat loss rate by the inelastic 
collisions. In the nearly elastic limit, the leading order contributions to the 
viscosity n and heat conductivity n are given by their classical values for for 
elastic hard disks: 

1 /T\ 2 / T ^ 



7? = — — and k = - — 

1 2d\TT ) d\TT 

The heat loss rate is the following: 

/ = 7r 1 / 2 (l-r 2 )dn 2 T 3 / 2 . 
The boundary conditions for the problem include the fixed temperature 

T(y = 0,t) = T (4) 

and zero velocity 

v(y = 0,t) = (5) 
at the base. Now, as the total number of grains TV is constant: 

poo 

/ n(y,t) dy — N/L x — const , (6) 
Jo 

the mass flux should vanish at infinity: 

nv = at y = oo. (7) 

Similarly, one should assume that the momentum flux at infinity is zero: 

_ 2 dv 
nl +nv — 77— = at y = oo, (8) 

as well as the energy flux: 

fv 2 ^\ dv dT n 

nv h 2T - 77 7J- k— =0 at y = oo. (9) 

\ 2 / dy dy 

We refer the reader to the book of Landau and Lifshitz [15] for a detailed 
discussion of the divergent forms of the hydrodynamic equations, and for the 
expressions for the momentum and energy flux density used in Eqs. (8) and 
(9). 

What are the characteristic length/time scales of the problem? The fixed 
base temperature To and gravity acceleration g define a characteristic macro- 
scopic length scale A = T /g. The characteristic number density of the 
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gas is therefore n = N/(XL X ). Correspondingly, there are three charac- 
teristic time scales in the problem. The fast time scale of the problem, 

1 /2 

Tf = T ' /g, is independent of the particle size. The heat conduction time 
scale is t = (-k\/ g) 1 / 2 (Nd/ L x ). The inelastic heat loss time scale 77 can be 
defined as t^ 1 = (1 - r 2 )dn T^ 2 . 

How do we know that the time scale 77 is indeed fast? Before we answer 
this question, let us first make a preliminary evaluation of the validity of hy- 
drodynamics. We should demand that the characteristic mean free path of 
the particles ~ (nod) -1 be much less than the macroscopic length scale A. 
This immediately follows L x /(Nd) -C 1. Notice that the quantity Nd/L x is 
of the order of the number of monolayers at rest (that is, when the system is 
not fluidized). Therefore, the validity of the hydrodynamics requires that the 
number of monolayers be much larger than 1 . Now we compare the character- 
istic time scales 77 and r introduced above. We see that the ratio 77 /r is of 
the order of L x /(Nd), which should be much less than unity if we want the 
hydrodynamics to be valid. 

To make a full use of hydrodynamics, we introduce scaled variables. Let 
us measure time in the units of the heat conduction time r, the coordinate in 
the units of A, the temperature in the units of To, the density in the units of 
no and the mean flow velocity in the units of A/r. In the scaled variables, the 
governing equations (l)-(3) become 



n ('-I- + v "pl + nT^ = fV/ 2 (»!) 2 + ££t»/» - 4AVT»/», (12) 



\dt dy J dy 2 \dy J ' 3 dy 2 
where 

L x (l-r 2 )V2 

£= 77^73 <1 and A=± — ' 

are the two governing scaled parameters of the problem. Parameter A is of 
the order of the square root of the ratio between the heat conduction time 
scale and inelastic heat loss time scale. 

The reduction of the four-dimensional space of governing parameters 
d, N, r and L x to the (e, A) - plane is an immediate benefit from using hydro- 
dynamics. Indeed, series of experiments or simulations with different d, N, r 
and L Xl but with the same e and A, should produce the same dynamics in 
the scaled variables. We will see in the following that, in the late stage of the 
dynamics (and in the steady state), the only relevant parameter is A. 

In the scaled variables, the boundary conditions (4), (5), (7), (8) and (9) 
become T = 1 and v = at y = 0; and 
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nv = 0, (13) 

nT + e 2 nv 2 ~ £ \ 1 / 2 ^=0 (14) 
2 dy 

and 

»»(^+r)-^„|-rV^o (15 , 

at y = oo. 

The simplest way to satisfy the boundary conditions (14) and (15) is to 
assume a zero heat flux at infinity, which implies dT/dy = at y = oo. As 
the total amount of material is finite, the density n should vanish at y = oo. 
Therefore, the nT term in Eq. (14) should also vanish there. Now, using 
what remains of Eq. (14) together with Eqs. (13), we see that Eq. (15) can 
be reduced to the simple condition dv/dy = at y = oo. Therefore, we 
are left with simple boundary conditions: T = 1 and v — at y = 0, and 
dT/dy = dv/dy = at y = oo [16]. A full formulation of the time-dependent 
problem requires that we specify the initial conditions n(y, t = 0),T(y, t = 0) 
and v(y,t = 0). We postpone this issue until later. 

Like in many other one-dimensional gasdynamic problems [17], it is conve- 
nient to make a transformation of variables from the Eulerian coordinate y to 
the Lagrangian mass coordinate m(y, t) = n(y' , t) dy' . The physical mean- 
ing of the Lagrangian mass coordinate is very simple: it is the mass content 
of the gas on the Eulerian interval (0, y) per unit length in the ^-direction. 
Notice that the infinite interval < y < oo in the Eulerian coordinates is 
mapped, by this coordinate transformation, into a finite interval < m < 1, 
at all times. Going over to the Lagrangian coordinate to, we can rewrite Eqs. 
(10)-(12) as follows: 



e 2 



d_ /1\ _ dv_ 

dt \n J dm ' 

dv d , „ N . e 2 d / i io dv 



(16) 



dt dm 2 dm V dm 



m-l + ^^-(T^n—\ (17) 



&T T dv _ g 2 rl/2 (9v\ 4 g ( 9 T 3/2\ 4yl 2 nT 3/2 /jgx 

dt dm 2 \ dm J 3 dm \ dm J 

The boundary conditions become T = 1 and v — at m = 0, and ndT/dm = 
ndv/dm — at to = 1. 

Notice that the transformation from the Eulerian to Lagrangian coordi- 
nate involves the density n(y, t) that is unknown a priori. We should not worry 
about it. After having determined the density n(m, t) in the Lagrangian co- 
ordinates, one can go back, at any time t, to the Eulerian coordinate by using 
the relation 

/ \ f m dm' 
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3 Steady state profiles and density inversion 

Let us briefly review the steady state solution of the problem [13, 14, 18]. To 
obtain the steady state profiles n s (m) and T s (m) we simply put d/dt = v = 
in Eqs. (16)-(18). Then Eq. (16) is obeyed automatically, while Eqs. (17) and 
(18) become 

4- (n s T s ) + 1 = (20) 
dm 

and 

'•' / d T^A - 3A 2 n s T^ 2 = 0, (21) 



dm \ s dm 

respectively. As we see, the steady state is described by a single scaled parame- 
ter A. Integrating Eq. (20) with respect to m and using the boundary condition 
at to = 1 (that becomes simply nT — 0), we obtain n s T s = 1 — to. Now we 
express the density through the temperature: n s (m) = (1 — m)/T s (m). Using 
this relation in Eq. (21), we arrive at the following linear equation [13, 14, 18] 
for uj s (m) = T^ 2 (m): 

(1 - m)u" -J s - A 2 {\ - m)uj s = 0, (22) 

where primes stand for the m-derivatives. This is a Bcsscl equation whose 
general solution is a linear combination of the modified Bessel functions of 
the first kind Iq [A(l — m)] and K [A(l — to)]. Demanding that the heat flux 
vanish at m — 1 and returning to the temperature variable, we obtain: 



Respectively, 

/x 1-to (l-m)In(A) .„„, 

To transform back to the Eulerian coordinate y, we should evaluate the 
integral 

dm' 



= / — 7 — a (25) 
Jo n s (m') 

numerically. Figures 2 and 3 show the steady state temperature and density 
profiles in the Lagrangian and Eulerian coordinates for A = \/l0. 

Equation (24) enables one to find the critical (minimum) value of A for 
the density inversion. At small A the density decreases monotonically with to 
(and therefore with y). The density inversion is born, at m = 0, at A = A c . 
Therefore, to find A c , we demand that the derivative n' s (m) vanish at to = 0. 
Using Eq. (24), we can reduce this condition to the equation Iq(A) = 2 A h(A), 
first obtained in Ref. [18]. Its solution yields A c ~ 1.06569 . . . . At A > A c the 
density reaches its maximum in the Lagrangian point to, = 1 — A c j A > 0. 
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Therefore, the position of the density peak in the Lagrangian coordinate grows 
with A monotonically. The respective (scaled) density maximum value is 

Hmax ~ AI$(A C ) • (26) 

Using Eq. (25), we can compute the scaled Eulcrian coordinate at which 
the maximum density maximum occurs: 

^T^T?^ (27) 



The dependence of on A is non-monotonic, see Fig. 1. Very close to the 
density inversion birth (at A > A c ) this dependence is linear: ~ A/A c — 1, 
while at A » A c one obtains y» ~ (2/1) Actually, the integral in Eq. (27) 
can be calculated analytically, in terms of a hypergeometric function. 



4 Formation of a density inversion: Low-Mach-Number 
flow 

Now let us return to the time-dependent problem. Can we exploit the presence 
of a small parameter e 2 in Eqs. (17) and (18) to simplify the equations? One 
can expect that the respective terms cannot be neglected during the rapid 
initial stage of the dynamics. The duration of this stage is of the order of the 
acoustic time or, in our scaled units, of the order of e. A comparison of different 
terms in the momentum equation (17) shows that the typical scaling of the 
velocity during this stage is e _1 3> 1. In the physical units, the flow velocity 
here is comparable to, or larger than, the speed of sound. Correspondingly, at 
early times the £ 2 -terms are comparable to, or larger than, the other terms 
and therefore should be kept. 

The situation changes at later times, determined by the inelastic heat 
losses and heat conduction. Here the characteristic Mach number of the flow 
becomes small, and we can neglect, in the leading order, the e 2 -terms in Eqs. 
(17) and (18). As the result, the momentum equation (17) is reduced to the 
hydrostatic condition (20), while in the energy equation (18) one can drop 
the viscous heating term. The continuity equation (16) remains the same as 
before. The boundary conditions at to = do not change, but those at to = 1 
are reduced to ndT/dm = 0. The boundary condition ndv/dm = should 
be dropped altogether. 

What is the physical meaning of the reduced model? A low-Mach-number 
flow proceeds on a time scale that is much longer than the time needed for 
establishing the (approximate) hydrostatic equilibrium, so the flow does not 
violate the hydrostatics. Similar reductions of nonlinear time-dependent gas- 
dynamic problems, that employ a separation of time/length scales, have been 
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successfully used in the context of the dynamics of an optically thin plasma 
that is cooled by radiation [19]. The presence of a heat loss term in the energy 
equation makes the latter problem similar to problems of granular dynamics. 

The hydrostatic condition immediately yields p(m,t) = n(m,t)T(m,t) = 
1 — m, like in the steady state solution! Using it together with the continuity 
equation (16), we can rewrite the energy equation as an evolution equation 
for u(m,t) = T 1 /^™, t): 

.„ , d 2 uj duj , 9 ,„ . . . 

u at = ^- m) a^-m- A ^- m)u - (28) 

This nonlinear parabolic equation should be solved with two boundary con- 
ditions: w(0,t) = 1 and (1 — m)du)/dm — > as m — > 1. Notice that the only 
relevant scaled parameter in the low-Mach-number problem is A. 
Once Eq. (28) is solved, we can find n(m,t) and v(m,t): 

, , d f m dm 1 
v(m,t) = — , , . . 29) 
v ; dtj n{m',t) v 1 

To transform back to the Eulerian coordinates we use Eq. (19). 

Equation (28) was solved numerically with a MATLAB PDE solver. Fig- 
ure 2 shows an example of numerical solution, at different times, for A = 
The initial condition for the temperature, T(m,t = 0) = 1, corresponds to 
a barometric density profile, which looks especially simple in the Lagrangian 
coordinate: n(m, t = 0) = 1 — m. The solution shows cooling of the gran- 
ular gas because of the inelastic heat losses (see Fig. 2 a). Simultaneously, 
a density peak develops at an intermediate height (Fig. 2 b). The material 
flux (Fig. 2 c) is everywhere negative: the material slowly falls from above 
until the steady-state with a zero mean velocity is reached. An important 
feature of the time-dependent solution is a transient temperature minimum 
that develops at intermediate times. It is caused by the fact that the inelastic 
heat losses are very small at large heights, where the gas density is small. 
Therefore, the gas cools down relatively fast at intermediate heights, which 
causes a temperature minimum there. Then a downward heat flux develops 
that removes the heat from above, finally creating a monotonically decreas- 
ing steady state temperature profile. One can see from Fig. 2 that the steady 
state temperature and density profiles coincide with the analytic solutions. 
The respective solutions in the Eulerian coordinates are shown in Fig. 3. 



5 Formation of a density inversion: Early times 

The low-Mach-number theory is invalid at early times, when the hydrostatics 
has not settled down yet, and the mean flow velocity can reach, or even ex- 
ceed, the local speed of sound. To investigate these early times, and to check 
the validity of the low-Mach-number theory at later times, we solved numeri- 
cally the full set of hydrodynamic equations (16)-(18) in the Lagrangian form. 
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The parameter e was taken to be \/5 • 10~ 2 ~ 0.02236. Equations (16) and 
(17) were solved with an explicit numerical scheme, based on the standard 
finite-difference method described in Ref. [20]. The viscous term in Eq. (17) 
was integrated explicitly, using a two-steps method, which is both stable and 
accurate. The energy equation (18) was treated implicitly, with a conservative 
scheme. 

We started with simple initial conditions: the constant temperature T(m, t = 
0) = 1, the barometric density profile n(m, t = 0) = 1 — m and a zero mean 
velocity. As expected, the late-time results show good agreement with the pre- 
dictions of the low-Mach-number theory. Figure 2 shows that, at latest times, 
the full-hydrodynamics solutions coincide with those obtained with the low- 
Mach-number theory. At an intermediate time one can see a local increase of 
the temperature at high altitudes, which is predicted by the full hydrodynam- 
ics, but missed by the low-Mach-number theory. The reason for this transient 
temperature maximum is the relatively large viscous heating that occurs at 
high altitudes because of a steep velocity gradient there. 

The early-time results of the full numerical solution, transformed back to 
the Eulerian coordinate, are shown in Fig. 4. One can see that the devel- 
opment of the density peak, and the formation of the transient temperature 
minimum, start already at early times. A more detailed analysis does show 
differences between the fast and slow theories. Indeed, a necessary condition 
for the validity of the low-Mach-number theory is the hydrostatic balance con- 
dition. This condition looks especially simple in the Lagrangian coordinate: 
p(m, t) = 1 — m. That is, when the hydrostatics sets in, the plots of the pres- 
sure versus m, at different times, should collapse into a single straight line. 
Figure 5 shows the pressure plots at different times. One can see that, at late 
times, the hydrostatic condition sets in as expected. At early times pressure 
oscillations are observed, and hydrostatics does not hold. Figure 6 shows the 
pressure history at the bottom plate m — 0, or y — 0. One can see that the 
pressure undergoes oscillations on the fast time scale of the system O(e). The 
oscillations rapidly decay with time and, starting from t ~ 0.4, the pressure 
stays almost constant which implies a low-Mach-number flow. 

6 Discussion 

We have employed granular hydrodynamics to describe a time-dependent 
granular flow that gives rise to a density inversion in a vibrofluidized granular 
system. We have predicted the formation of a transient temperature minimum 
at intermediate times. The temperature minimum ultimately disappears as 
the system approaches the steady state with a density inversion. We have also 
predicted pressure oscillations at early times, when the hydrostatic condition 
is not yet satisfied. 

We hope that the predicted phenomena will be observed in molecular 
dynamics (MD) simulations and in experiments with vibrofluidized granular 
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beds. The MD simulations should start from the barometric density profile, 
corresponding to the temperature of the bottom plate. A natural way to initi- 
ate the initial state is to first run the simulation with elastic particle collisions. 
After the barometric profile sets in, the collision rule should be switched into 
inelastic regime. 

The hydrodynamic theory considered in this work is expected to be quan- 
titatively accurate for nearly elastic collisions, 1 — r 2 <C 1, and in the dilute 
limit. The latter condition demands that the maximum number density of the 
particles be less than, say, 0.2 n c , where n c = 2/(v / 3cP) is the hexagonal close 
packing density. The maximum density peak is achieved at the steady state. 
Therefore, using Eq. (26) and going back to the "physical" density, we can 
write the dilute regime criterion as 



A c V3Nd 2 I$(A) 



< 0.2 . (30) 



I'(A C ) 2\L X A 
An alternative form of the same criterion is 

^£7 < 0.65. (3D 

e A A y ' 

We expect that the theory is still valid qualitatively when this criterion is 
not satisfied, but the maximum density is still less than about 0.5 n c . To 
achieve a better accuracy for such moderately dense flows, one can use the 
constitutive relations of Jenkins and Richman [11] that account for excluded 
volume effects. 
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Fig. 1. The density peak height in the (scaled) Eulerian coordinate versus A. 
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Fig. 2. The temperature (a), number density (b) and mass flux nv (c) versus the La- 
grangian mass coordinate m at late times, as computed from the low-Mach-number 
model (symbols) and full hydrodynamics (lines) . The scaled parameters are A = VlO 
and e = \/Z ■ 1CT 2 . The scaled times are t = 0.1 (solid circle, dash-dot-dot line), 0.3 
(solid square, dashed line), 0.5 (empty circle, dotted line) and 0.85 (empty square, 
solid line). The analytic steady state solutions T a (m) and n s (m) are indistinguish- 
able from the respective solid lines. 
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Fig. 3. The temperature (a), number density (b) and mass flux nv (c) versus the 
Eulerian coordinate y at late times, as computed from the low-Mach-number model 
(symbols) and lull hydrodynamics (lines). The scaled parameters are A = VlO and 
£ — V5 ■ 1CP 2 . The figures on the left and on the right show the same profiles on 
different scales ol height. The scaled times are t = 0.1 (solid circle, dash-dot-dot 
line), 0.3 (solid square, dashed line), 0.5 (empty circle, dotted line) and 0.85 (empty 
square, solid line). The steady state solutions T s (y) and n B {y) are indistinguishable 
from the respective solid lines. 
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Fig. 4. The temperature (a), number density (b) and mass flux nv (c) versus the 
Eulerian coordinate y at early times, as computed from the full hydrodynamics. The 
scaled parameters are A = \/T0 and e = y/E- 10~ 2 . The scaled times are t = 0.01 (the 
dotted line), 0.04 (the solid line), 0.07 (the dash-dot-dot line) and 0.1 (the dashed 
line) . 
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Fig. 5. The pressure p = nT versus the Lagrangian mass coordinate m at different 
times, as computed from the full hydrodynamics. The scaled parameters are A — 
y^W and s = VE ■ 1CT 2 . The scaled times are t = 0.01 (the dashed line), 0.04 (the 
dash-dot-dot line), 0.5 (the solid line) and 0.85 (the empty circles). 
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Fig. 6. The pressure p — nT at the bottom plate y = 0, versus time, as computed 
from the full hydrodynamics. The scaled parameters are A = \/l0 and e = \/E- 1CP 2 . 



